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ABSTRACT 

Prior to the launch of the Swift mission several X-ray line detections were 
reported in Gamma Ray Burst afterglow spectra. To date, these pre-Swift era 
results have not been conclusively confirmed. The most contentious issue in this 
area is the choice of statistical method used to evaluate the significance of these 
features. In this paper we compare three different methods already extant in 
the literature for assessing the significance of possible line features and discuss 
their relative advantages and disadvantages. The methods are demonstrated by 
application to observations of 40 bursts from the archive of Swift XRT at early 
times (< few ks post trigger in the rest frame of the burst). Based on this 
thorough analysis we found no strong evidence for emission lines. For each of the 
three methods we have determined detection limits for emission line strengths in 
bursts with spectral parameters typical of the Swift-eia sample. We also discuss 
the effects of the current calibration status on emission line detection. 

Subject headings: gamma rays: bursts — gamma rays: observations — line: 
identification — methods: data analysis — methods: statistical 
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Introduction 



It is widely accepted that the spectra of the X-ray afterglow of Gamma-Ray Bursts 
(GRBs) are domin ated by non -thermal emission, the leading candidate for which is syn- 
chrotron emission (jPiraru |2005| and references ther ein), though alt ernat e emission processes 



have also been suggested such as self-Compton (IWaxman 



199 



3 and 



1999r) or inverse Compton scattering o f external light (IBrainerd et al. 
Shaviv fc Darlll995l : llazzati et al1l2004l ). 



Ghisellini fc Celotti 



1994 



Shemi 



1994 



Up to the present time the X-ray spectra of Sw ift afterglows are generally well described 
by an absorbed power law (for counter-examples see iButlerl 120071 ). typically absorbed by ma- 



terial with a colum n density in exce s s of th e well measured Galactic values (jCampana et al. 



20064). Table 2 of ICampana et al.l (l2006d ) shows that, of the 17 bursts analyzed, 14 have 
observed Nn values greater than the measured Galactic column density, whilst the remaining 
three have observed A^h values that are consistent, within limits, with the measured values. 

In the past it has been proposed that there are other spectral featur es, with varying 



levels of significance, in addition to the 


3asic absorbed power law spectrum ( 


Piro et al. 


1999; 


Yoshida et al. 


1999; 


Amati et al.l 


2000 


Antonelli et al.l 


2000 




Piro et al. 2000 




Reeves et al. 


2002; 


Watson et al. 


2002;IWatson et al. 


2003 and 


Frontera et al. 


2004 


). Most 


are attributed 



to Fe Ka emission lines or the radiative recombination continuum of the same element. Some 
have been attributed to the Ka lines of Ni, Co or of lighter elements such as Si, S, Ar and 
Ca. In tw o cases there has been a report of a trans ient absorption feature also corresponding 
to Fe jAmati et allboool : IProntera et al1l2004l ). 



The models for the production of such emission features are divided into transmission 
and reflection models, though the large equivalent widths (~ few keV) infer red from the ob- 
served X-ray features favor models i n whi ch th e line is produced b y reflection (IRees fc Meszaros 



2000l : iBallantyne fc Ramirez- Ruia 1200 ll and IVietri et al.l 1200 ll ). Proposed models have to 



overcome two constraints; the size problem and the kinematic problem. Observing a line at 
a time tobs after the burst implies that the emitting material must be within a distance of 
~ ctobs/(l + z) from the central engine, thus implying that the region must be compact if 
a line, or lines, are observed at early times. Additionally the emitting region must contain 
~ O.IMq of Fe (in the case of Fe features) whilst still being optically thin to electron 
scatteri ng, in order that Comptonization does not broaden the line beyond the observed 
widths ( IVietri et al.ll200ll ). If the line width is interpreted as being due to the velocity of the 
supernova remnant, the observed limit on this width implies an age limit on the remnant of 
~ 10 — 20 days. However, at this time, Co nuclei outnumber both Ni and Fe nuclei; thus the 
emission line would be due to Co at an energy of 7.5/(1 + z) keV, which is the kinematic 
problem. 
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Various geometries have been suggested for the reflection models, which rely on either a 
precursor or simultaneous supernova (SN) event. If a SN occurs several tens of days before 
the GRB this solves both the size and kinematic problems. In these cases the radiation 
from the GRB jets can either illuminate the inner face of the SN shell remnant or the inner 
faces of wide funnels that they excavate through young plerionic remnants. However, these 
models have been questioned follo wing the simultaneous GRB-SN assoc iation indicated by 
GRB 980425 jGalama et al.lllQQsh and then confirmed by GRB 030329 jHiorth et al.l l2003l: 
Stanek et all l2003h and~GRB 060218 JCampana et al.l boOGbl : IPian et al.l bood l In this 
case the most likely scenario for emission line production occurs if the progenitor ejects a 
large amount of matter, at subrelativistic speeds, along its equator. The halo of material 
surrounding massive stars, ejected by their strong stellar winds towards the end of their main 
sequence lifetime, scatters a fraction of the photons from the prompt and afte rglow phase 



back into the equatorial material, which then produces X-ray line emission (jVietri et al. 
200lh . 



Verifying the presence of such spectral features is of critical importance as they will 
allow us to probe the circumburst environment of the GRB as well as gaining an indirect 
indication of the possible structure and behavior of the central engine. 

The statistical significance of the 1999-2003 reported featur es is low (usuall y 2 — 3a), 
only two detec tions have a significa nce >4cr (GRB 991216: 4.7a, iPiro et al.ll2000l : and GRB 
030227: 4.4cr, IWatson et al.ll2003l ). Though later detections were made with much more 
sensitive instruments than the early ones, all detections have remained at this low significance 
level, as a result they remain the subject of much debate. Arguably the most controversial 
issue in the discussion of line detections is the choice of statistical method employed to gauge 
their significance. At least four methods have been proposed and used in the GRB literature: 



1. The likelihood ratio test (LRT) and related F-test 

2. Bayes factors 

3. Bayesian posterior predictive probability 

4. Monte Carlo test for peaks in data after 'matched filter' smoothing 



ication of these me t hods to GRB X -r ay spectra can be found in: 



Freeman et al.l ( 


1999).^ 


"t^oshida e 


t al. 


(1999). |] 


^iro et al. 


(1999 


). 


Protassov et al. 


J2OO2). 


Rutledee & Sako 


(2003 




Tavecchio et a] 


. ( 


2004 


), 


Butler et al. 


( 


2005). J 


5ako et al. 


( 


2005). 


Butlei 


( 


2007 


) and 



references therein. 
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In all the applications cited above, an underlying continuum was assumed, usually 
in the form of an absorbed p ower law (e.g. using the Wisconsin absorption model, see 
Morrison fc McCammonI Il983l ). The detection of a line then amounts to a comparison of 
two models: Mq, the simple "continuum" model, and Mi, the more complex "continuum 
+ line" model. The strength, location and width of the emission line may be restricted or 
allowed to be free parameters. 



As discussed in depth by iProtassov et al.l (120021 ) and iFreeman et al.l (119991 ) , there are 
strong theoretical reasons why the LRT is not suitable for assessing the significance of emis- 
sions lines, despite its popularity in the literature. We will not repeat those arguments here. 
It is the purpose of the present paper to compare the relative merits of the remaining three 
methods, in terms of their computational efficiency, robustness and sensitivity limits by ap- 
plying all three methods to X-ray spectra from the Swift archive. This is a particularly rich 
archive because of t he combination of the rapid slew respon se of the 6*^;^^ GRB mission 
(IGehrels et al.l 120041 ) and the powerful X-ray Telescope (XRT; iBurrows et al.ll2005al ). 

The remainder of this paper is organized as follows. In §2]we provide details of the sample 
selection criteria and basic data reduction. The theoretical basis and practical applications 
of the three statistical methods under investigation are described in detail in ^ In §1] we 
apply all of the methods to PKS 0745-19, a known line emitting source, to demonstrate the 
expected outcome when a line is present, and [|5] discusses a simulation study to assess the line 
detection limits of each method for typical Swift XRT data. In ^Hlwe discuss the results from 
the Swift archival GRB afterglow data, highlighting several GRBs with potential additional 
spectral components, ^is dedicated to a discussion of our results and their comparison to 
other recent line detections in the literature. Finally, ^ presents our conclusions. 



Data Reduction 



This paper reports on the analysis of Windowed Timing (WT) mode data from GRB 
050128 to GRB 060510B, covering a total of 153 bursts, 40 of which contained sufficient 
WT mode data for our analysis methods. WT mode data was chosen primarily because the 
time interval covered by these observations, typically T-l-0 s to T-|- 500 s (though for bright 
bursts this may extend up to T-|-few ks) in the rest frame of the burst, is rarely explored. 
Prior observations of the 1999-2003 bur sts typically sta r t at 2 0-|- hours after the trigger 
in the observer's reference frame, though lAntonelli et al.l (120001) report o n an emission line 
detection at T+~ 12 hours; lAmati et al.l (120001 ) and lFrontera et al.l (120041 ) report absorption 
line features in very early time data (T-|-<20 s and T-|-<300 s respectively) from the WFC 
of BeppoSAX . Additionally WT mode data is only taken whilst the GRB afterglow is bright. 
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All of the methods discussed can easily be extended to Photon Counting (PC) mode data. 
We acknowledge that the current theoretical models for line emission indicate that lines 
could occur at times not covered by WT mode data, however, the same models do not rule 
out this time period either. 

All data have been obtained from the UK Swift archiv^ ( Tyler et al. 2006 ) and pro- 
cessed through xrtpipeline vO.lO.jj^ using version 008 calibration files and correcting for the 
WT mode gain offset (if p resent). Version 008 o f the CALDB is a marked improvement 
over the previous release ( Campana et al. 2006a however, it may be the case that low 
energy calibration features have still not been optimally corrected. We did not apply any 
systematic correction factor to the errors of our spectra because the recommended factor is 
very much smaller than the statistical errors in our spectra. Grade — 2 data, using ex- 
traction regions of 20x3 pixels, for both source and background regions have been used. At 



count rates below 100 counts s ^ WT mode data does not suffer from pile-up (IRomano et al. 



20061 ). however, s ome of the time interv als considered contain sufficient flux to cause pile-up 
effects. Following iRomano et al.l (120061 ) we have excluded central regions when necessary as 
detailed in their Appendix A, splitting the 20x3 pixel region into two 10x3 pixel regions 
placed either side of the central exclusion region. 

All spectral fitting and simulations have been carried out using XSPEC version 12.2. lab 
or higher with background subtracted spectra binned to >20 counts bin~^. This binning 
permits the use of the minimization as a Maximum Likelihood method. Data for each 
GRB being considered have been time-sliced with the following criteria in mind: 



1. Each spectrum must contain 800 — 1600 background subtracted counts, 
compromise between good time resolution and spectral quality. 



This is a 



2. If one or more flares are present in the data, wherever possible (not violating condition 
1), separate spectra were extracted for the risin g and falling sections of the flare, s ince 
spectral evolution is likely to occur at this time (IBurrows et al.ll2005bl . IZhang||2007l and 
references therein). 

3. If data are affected by pile- up then these time periods were extracted separately from 
the non piled-up data. 



^http://www. swift. ac.uk/swift_livc/obscatpage.php 
^Release date 2006-03-16 

3http://hcasarc.gsfc.nasa.gov/docs/heasarc/caldb/swift/docs/xrt/SWIFT-XRT-CALDB-09.pdf 
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The range of 800 — 1600 counts was chosen to ensure good time resolution while main- 
taining sufficient counts to obtain a reasonable spectrum, with ~ 40 — 80 spectral bins (each 
with ~ 20 counts) over the useful bandpass. The data considered here were taken during the 
early, bright phases of the afterglow evolution (typically T+0 s to T+500 s), during which 
the X-ray flux and (possibly) spectrum are often highly variable. Time resolution is there- 
fore important to reduce the effects of flux/spectral variation on the modelling of individual 
spectra. Furthermore, previous claims of emissions lines have often reported the features as 
transient and so time resolution may be important for detecting lines. 



3. Analysis Methods 

As noted in ^ several different methods have been used in the past to assess the sig- 
nificance of line detections in the X-ray spectra of GRBs. The three methods that are the 
subject of the present paper are discussed individually in the following subsections. The 
reader interested only in the application of these methods may wish to skip to § HI 



3.1. Bayes Factors 

The goal of scientific inference is to draw conclusions about the plausibility of some 
hypothesis or model, M, based on the available data D = {xi,X2, ■ ■ ■ ,xn}, given the back- 
ground information I (such as the detector calibration, statistical distribution of the data, 
etc.). However, when presented with data it is usually not possible to compute this directly. 
What can be calculated directly in many cases is the sampling distribution for data assum- 
ing the model to be true, p{D\M,I). This is usually called the likelihood when considered 
as a function of M for fixed D. Statements about data conditional on the model may be 
related to statements about the model conditional on the data by Bayes' Theorenj^. In its 
usual form Bayes theorem relates the likelihood to the posterior probability of the model M 
conditional on data D (and any relevant background information /), written p{M\T), I): 

The term p{M\I) is the prior probability of the model M and describes our knowledge (or 
ignorance) of the model prior to consideration of the data (often called simply the 'prior'). 
The term p(D\I) is effectively a normalization term and is known as the prior predictive 



^For general references relating to Bayesian analysis see 



http://www.astro.cornell.edu/staff/loredo/bayes/ 
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probability (it describes the probability with which one would predict the data given only 
prior i n formation abo ut the model). F o r a more gen e ral discussion o f Baye s theorem s ee iLee 
fll989h . Ihoredol Jl99oh . Ihoredol Jl992h . Isivial Jl996h . ICelman et all Jl995h. ICregorvl J2005h 



and for discussion in the co ntext of GRB line searches see lFreeman et al.l (119991 . their §3.1.2) 
and iProtassov et al.l (120021 ). In the rest of this paper we drop the explicit conditioning on 
background information /, but it is taken a s accepted that "no probability judgements can 
be made in a vacuum" (IGelman et al.lll995l ). 



One simple way to represent the poster ior probabilitie s for two alternative models is in 
terms of their ratio, the posterior odds (see lGregoryll2005l § 3.5). This eliminates the p(D) 
term (which has no dependence on M). If we define two competing models, such as one with 
a line (Mi) and one without (Mq), we may compute the posterior odds: 



O 



p(Mi|D) p(Mi)p(D|Mi) p{Mi 



10 



p(Mo|D) p(Mo)p(D|Mo) p(Mo 



-B 



10- 



(2) 



High odds indicate good evidence for the existence of a line in the spectrum. The first term 
on the right hand side is the ratio of the p riors, the second t erm i s the ratio of the likelihoods 
and is often called the Bayes factor (see iKass fc Rafteryl (119951 ) for a detailed review). In 
the present context we have no strong theoretical grounds to prefer one or the other model 
(line or no line) and so assign equal prior probabilities to our two models. Thus the ratio of 
the priors in equation [2] is set to unity and the posterior odds are equal to the Bayes factor. 
In the following we use the terms posterior odds, odds and Bayes factors interchangeably. 

The likelihood functions in equation [2] are functions of Mj only. If the models contain no 
free parameters (i.e. are completely specified) then equation [2] can be used directly. However, 
if the model does contain free parameters, the likelihood will be a function of the parameter 
values. In the present context, where the particular values of the parameters are not the 
subject of the investigation, the parameters are referred to as nuisance parameters. In order 
to remove the dependence on these nuisance parameters the likelihood function must be 
written as a function of the N parameters (denoted 6 = {6i, 62, ■ ■ ■ , On}) and integrated, or 
marginalized, over the prior probability density function (PDF) for the parameters: 



p(D\M)= p(D,0\M)de = p(D\e,M)p{e\M)de 



(3) 



The marginal likelihood is obtained by integrating over all parameter values the joint PDF 
for the data and the parameters. This joint PDF may be separated into the product of 
two terms using the rules of probability theory: p(D\0, M) is the likelihood function of the 
data as a function of the model and its parameters, and p{6\M) is the prior for the model 
parameters. Once these are assigned one can compute the necessary likelihood (a function 
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of M alone) by integration. The Bayes factor for model Mi (with parameters di) against 
model Mq (with parameters ^o) iiiay now be written: 

/p(D|gi,Mi)p(0i|Mi)rf^i 
'° /p(D|0o,MoM0o|Mo)rf0o' ^ ^ 

The issues of how the relevant likelihoods and priors are assigned, and the integrals computed, 
are discussed below. 



3.1.1. Application to high count X-ray spectra 

In the limit of a large number of counts per spectral bin, the Poisson distribution of 
counts in each bin will converge to the Gaussian dist ribution, and in t his case equation [3] 



can be written in terms of the familiar x fit statistic (lEadie et al.lll97ll ): 

1 ^ ^ ( — (RW"^ 

L = \n[piB\e, M)] = ~Y1 1^[2^^'] - E 2a' = ~ ^'/^ 

1=1 1=1 « 

where cxj is the error on the data (e.g. counts) measured in the ith channel, and fi{6)i is the 
predicted (e.g. model counts) in the channel based on the model with parameter values 0. 
The last equality can be made since the term l^[27ro"j^] is constant given data with errors o"j. 
This is why, in the high count limit, finding the parameter values at which is minimized 
is equivalent to finding the Maximum Likelihood Estimates (MLE) of the parameters, 6. 



3.1.2. Approximating the posterior 



In general the integrals of equation [3] and H] must be computed nuraerical ly, using for 
example Mark ov Chain Monte Carlo (MCMC) methods (iGelman et al.l Il995l chapter 11; 
Gregoryl |2005| chapter 12), which is computationally demanding. However, maximum like- 
lihood theory says that the MLE will become more Gaussian and of smaller variance as 
th e sample size (number of counts) increases, even if the model is non-linear (chapter 11 
of iGregoryl |2005| ) . Therefore, with sufficient counts the likelihood function will approach a 
multidimensional Gaussian with a peak at the MLE location 6, i.e. the location of the best 
fit (minimum x^) iii parameter space. Furthermore, if the prior function is relatively fiat 
around the peak of the Gaussian likelihood we may approximate the prior term in equation [3] 
by a constant, namely its value at the best fit location, p{6\M). Putting this together means 
we may approximate the posterior as a Gaussian - often called the Laplace approximation 
- which greatly simplifies the integrations in equations [3] and HI since a multidimensional 
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Gaussian may be evaluated analytically, once its peak locations and covariances are known, 
which avoids the need for computationally expensive numerical integration. 



The integral of an unnormalized multidimensional Gaussian is (27r)^/^A/det [cr^] times 
the peak value, where cr^ is the covariance matrij|§ evaluated at the peak (best fit location) 
and N is the number of parameters. We may now re- write equation [3] as: 



(6) 



which involves the prior density only at the mode of the p{6i\Mi) likelihood (i.e. the density 
at the MLE posi tion) for each model. This can be substituted into equation [2] to give the 
Bayes factor (see lGregoryll2005l chapters 10-11): 



Bio = 



e.p(-A,V2)(2.)-'^«^ 

Vdetfcr^] p(6»o|Mo) 



(7) 



Xi — Xo- This can be calculated, using the appropriate 



where AA^ = Ni - Nq and Ax^ 
values of and the covariance matrix cr^ evaluated at the best fit location for each model 



once we assign prior densities to each parameter. See iGregoryl (120051 ). chapters 10-11 for a 
discussion of essentially the same method. 



3.1.3. Validity of the Laplace approximation 

There are a number of ways to check the validity of this assumption. One is to inspect 
the shape of the x^{0) surface, which is related to the likelihood surface by L ~ — x^/2. A 
Gaussian likelihood is equivalent to a paraboloidal log-likelihood or If contours 

of Ax^ appear paraboloidal around the mimimum, in one and two dimensions, for each 
parameter or pair of parameters, the likelihood surface must be approximately Gaussian (the 
conditional and marginal distributions of a Gaussian are also Gaussian, so slices through the 
X^ space should also be paraboloidal). This was generally true for the continuum model for 
the XRT data. 



^ The covariance matrix is the square, symmetric matrix comprising the covariances of parameters 9i and 
9j as element afy By symmetry af^ = cr|j. The diagonal elements are the variances of the parameters. The 
covariance matrix may be estimated as minus the inverse the Hessian matrix listing all the second derivatives 
of the log likehhood function [VVL]y = d'^L/dO^dOj. Given that L = log[p(D|0, M)] = constant - (in 
the limit of many counts per channel) the covariance matrix may be estimated using [cr^Jij — 2[(VVx^)~^]ij- 
The second derivatives of the X^(^) function can be evaluated numerically. 
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As a further test of the Gaussian approximation we compared the posterior calculated 
using the L aplace approximation with the posterior calculated using the MCMC algorithm 
discussed in Ivan Dyk et al.l (1200 ll ). The MCMC method does not use an analytical approx- 
imation for the posterior, and therefore is a more general method, but is computationally 
demanding. Figure [1] illustrates the two posterior distributions calculated for the specific 
case of a spectrum from GRB 060124. The Gaussian data were computed from 10^ random 
draws from a multidimensional Gaussian with a covariance matrix evaluated as the minimum 
location using XSPE C. The non-Gaussian d ata were generated from 10^ draws generatecl§ 
by the MCMC routine of Ivan Dyk et al.l (1200 ll ). It is clear that the two distributions are not 
identical but are very similar both in terms of size and shape. In the present context it is 
important that the "credible regions" occupy similar volumes of parameter space. 

The above analyses demonstrate the Gaussian approximation is reasonable for the pos- 
terior of the simple "continuum" model Mq, which is the denominator of equation HI The 
same will be true of the mo re complex "co ntinuum + line" model Mi when the line is well 
detected (see section 11.3 of lGregoryll2005l ). When the line is weakly detected the posterior 
will be close to the boundary of the parameter space, in which case the Gaussian approx- 
imation will not be so accurate. Indeed, when the MLE of the line normalization is close 
to the boundary the likelihood (and therefore posterior) enclosed in the allowed region of 
parameter space will be smaller than that given by the Laplace approximation, which as- 
sumes the Gaussian function extends to infinity in all directions. This will also happen, for 
example, when the best-fitting line energy is near the limit of the allowed energy range. In 
such cases there will be a tendency to overestimate the Bayes factor (i.e. favor Mi). But 
when the line is weak there may be multiple peaks in the likelihood (and posterior) which are 
not accounted for explicitly in the Laplace approximation. We therefore treat the calculated 
Bayes factors only as a rough guide to the presence of a spectral line. 



3.1.4. Assigning Priors 

Bayes factors are sensitive to the choice of prior density. As stated above, using the 
Laplace approximation the resulting Bayes factors are sensitive to the prior densities only 
at the MLE parameter values, but we must exercise care in assigning prior density functions 
in order that these values are reasonable. Fortunately, the prior densities for all parameters 



^ FoUowing lvan Dyk et al.l (|2001r ) we generated five seperate chains, starting from different, 'overdispersed' 
positions within the parameter space (all outside the 99% confidence region calculated using Ax^) and used 
the R^/^ statistic to assess their convergence. We collected data from the chains only after R^/^ < 1.01. 
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that are common to Mq and Mi (such as photon index and normahzation) are the same for 
Mq and Ml, and therefore cancel out in the ratio. For the other parameters we have no 
cogent information except for their allowed range s . In such cases we should use the "least 



informative" prior densities (see e.g. iLoreddllQQOl : ISivialll996l : lGregoryll2005l and references 
therein for further discussion). 

There are wide ranges of possible line energies and redshifts and so the line energy, 
Eiine is only constrained to he within the useful XRT bandpass, typically 0.3 — 10 keV. We 
therefore assigned a uniform prior density p(i?iine| Mi) = 1/ [i^max — -E'min] over this range. For 
most spectral fits the line width W was initially held fixed at a value below the instrumental 
resolutior0, and later allowed as a free parameter. For those models in which the width of 
the line was a free parameter, the width was assigned a uniform prior over the allowed range 
(usually 0.0 - 0.7 keV): p{W\M,) = 1/[W^^^ - W^^^]. 

In order to test the dependence of the results to the prior densities, two non-informative 
pr ior assign r nents were made for the line strengt h (no r maliza tion) , A, following the discussion 



m 



Gregory! (120051 ). Firstly, following §4.2 of ISivial (119961 ) the line strength was assigned 
a uniform prior between zero and some upper limit Amax- Previous reports of ernission 



lines have estimate d the line flux A to be as lit tle as a few percent (jReeves et al.l 12002 
Watson et all boosh or as much as ~ 40 - 80% JVoshida et all Il999l : IPiro et all boooh of 
the total flux. We conservatively take Amax to be the total flux of the spectrum over the 
evaluated bandpass (i.e. our constraint is that the liiie flux is between 0-100% the source 
flux) . However, there are strong arguments (iLoredd Il990l : iGelman et al.l Il995l : [Gregory 
20051 ) that such as 'scale' parameter should be given a Jeffreys prior, p{A\Mi) ~ 1/A, which 
corresponds to a constant density in log(yl). Formally this is an improper prior (cannot be 
normalized such that its integral is unity), but one can apply reasonable upper and lower 
bounds in order to form a proper prior density. Following equation 3.38 of [Gregory] (120051 ) 
we used p{A\Mi) = l/y41n[AinaxMmin]- In the present context AmaxMmin = 800, since a 
reasonable lower limit to the X-ray counts from a line is one count, and a reasonable upper 
limit is 800, the total number of counts in the spectrum. This yields p{A\Mi) = 1/6. 68A as a 
normalised Jeffreys prior. The prior density is therefore higher for weaker lines in the Jeffreys 
case compared to the uniform case at values (i.e. over 1.25 x 10~^y4jnax < ^ < OAbA^i^^^), 
and lower for stronger lines. 



The ratio of the prior densities at the modes of the two likelihood functions is then 



^cr — 59 eV (at 5.895 keV) at launch (A. Beardmore, private communication) 
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simply 

piO,\M, 



p{Oq\Mo 



J9(EHne|Mi)p(H^|Mi)j9(A|Mi) 
1 

max -^min ][W^ax-W^min]Ap) 



max 



in the ranges i^une G [-Emin, -^max], G [W^mm, W^max] and zero elsewhere. Here, Ap = A 
in the uniform case or = Q.Q8A in the Jeffreys case. We have used both uniform and 
Jeffrey's priors in the analysis discussed below (see section IB. 111) . 



3.2. Posterior Predictive p-values (ppp) 

The use of posterior predic tive p- values (ppp) wa s advocated, and demonstrated by 



application to GRB spectra, by iProtassov et al.l (|2002| . see §4.1 for a description of their 



method and §5 for its application to GRB 970508). Like Bayes factors this method is 
grounded in Bayesian probability theory. 

One uses the posterior density, p(0|D), for the model parameters conditional on the 
data - which defines our state of knowledge about the parameters given the data and the 
available prior information - to determine the posterior predictive distribution - which is 
the distribution of possible future data predicted based on the observed data. ('Predictive' 
because it predicts possible future datasets and 'posterior' because the parameters are drawn 
from the posterior density of the parameters.) The posterior predictive distribution is: 

p(D^''^|D) = J p(D^™,6>|D)rf6> = J p{'D'''^\0)p{0\B)dO (9) 

where D*^'™ are the possible future datasets (simulations). In practice the posterior density 
is used to generate a set of random parameter values Of^ [i = 1,2,...) and each of these 
is used to simulate a random dataset Df™. The set of simulated data from all the possible 
random parameters defines the posterior predictive distribution for simulated data. This in 
turn can be used to define the posterior predictive distribution for some test statistic T(D) 
(which is a function of the data): 

p[T(D™)|D] = j p[T{B'''^)\e]p[0\B]d0 (10) 

(compare with equation [9]). The posterior predictive p- value (ppp) is the fraction of this 
distribution for which T(D™) > T(D), i.e. the area of the tail of the distribution with 
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values of the test statistic more extreme than the value from the observed data. 



/t(d) 



where the integration is taken over the posterior predictive distribution of D^'™. As such the 
ppp value is a Bayesian analogue of the p- value of null hypothesis tes t s fam iliar from classical 



statist ics (e.g. the F or tests). See chapter 6 of iGelman et al.l (119951) o r lGelman et al. 



( 1l996l ) for a general discussion of the ppp method, and lProtassov et al.l (120021 ) for application 
to GRB data. 

Using the posterior predictive distribution from equation |9] one can produce a large 
number of random simulated datasets to be used in a Monte Carlo scheme to calculate the 
integral of equation [11] numerically. The steps for a Monte Carlo method for computing the 
posterior predictive distribution to calibrate the test statistic T are as follows: 



1. Compute the value of the test statistic for the observed data, T(D) 

2. Randomly draw N sets of Mq model parameter values 6i for i = 1,2, ... ,N according 
to the appropriate posterior distribution p{6\D) 

3. For each of i = 1, 2, . . . , simulate a dataset Df™ using the randomly drawn parameter 
values 0i. This accounts for uncertainties in the parameter values. 

4. For each of the simulated datasets compute the test statistic T(Df™). This is the 
posterior predictive distribution of the test statistic given the observed data D. 

5. Compute the posterior predictive p- value as the fraction of simulated datasets that 
gave a test statistic more extreme than that for the observed data: 

1 ^ 

p = j;^J2^[nr>r)-Tm (12) 

i=l 

where G is the Heaviside step function which simply counts instances where T(Df™) > 
T(D). 



The number of simulations, A^, must be large to ensure a good approximation to the integral 
of equation [TT] (^ which is a multiple int e gral, being itself the integral of the function computed 
by equation [To]). See jProtassov et al.l (120021 ) for more detailed discussion. 
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3.2.1. Application to GRB X-ray spectra 

As discussed above we may approximate the posterior density for the parameters, using 
a multidimensional Gaussian centered on the MLE values and with a shape defined by the 
covariance matrix evaluated at the peak (<t^). The randomised parameter values needed for 
step 2 above may then be generated with the Cholesky method. 

For the purposes of the present paper we use as the test statistic the change in the 
fit statisti(Jf| between the two mode ls, Mq and Mi. This is equivalent to the formulation 



discussed in iProtassov et al.l (120021 ) . The observed data were fitted with the model Mq 
and the covariance matrix evaluated at the best-fit point used to construct the multivariate 
Gaussian distribution from which parameter values were randomly drawii^. For each set of 
model Mq parameter values a spectrum was simulated with the appropriate response matrix 
and exposure time, with counts in each channel drawn from a Poisson distribution, and 
binned in the same manner as the observed data. 

In order to calculate the test statistic for each simulation, T(Df™) it was necessary to 
fit each simulated dataset with the two competing models Mq and Mi, for each one find 
the best-fitting parameters, and compute Ax^. This necessarily involves a computationally 
expensive multi-dimensional parameter estimation for each of the simulations. We use 
as standard = 10"^ simulations which yields a p-value accurate to four decimal places at 
very highest and lowest p-values (there is an uncertainty on the ppp value from the finite 
number of simulations which is roughly a/p(1 — p)/N from the binomial distribution). This 
is acceptable for determining p-values as low as p ~ 10~^, i.e. 99.99% 'significance'. 

As a further test of the validity of the Gaussian assumption for the posterior (see also 
sections I3.1.21|3?T73|) we have compared results with and without this assumption. In partic- 
ular, we calculated the ppp-\a\\xe for a spectrum of GRB 060124 using Gaus sian parameter 



values and also using values generated by the MCMC method discussed by Ivan Dyk et al. 



( I2OOII ). The two results were reasonably close {p = 0.050 ± 0.007 from the Gaussian simula- 



tions and p = 0.073 ± 0.008 from the MCMC, based on 10^ simulations). This confirms the 



^ The Ax^ statistic is familiar to most X-ray astronomers and was used in the Bayes factors method 
above. Here we note that it is equivalent to the likelihood ratio test (LRT) statistic, since using equation [5] 



we have Ax = X(o) ~ X(i) = — 21nA, where A = pCDIOq, Mq) /p(T>\9i, Mi) is the ratio of the hkelihood 
maxima of the two models. Under the assumptions for which the LRT is valid this should be distributed as 
X^ with degrees of freedom equal to the number of additional free parameters in model Mi compared to Mq. 
The reason for chosing the LRT over relate d , stati stics, such as the F-test, is that LRT is more powerful. See 
Freeman et aD (|l999h and IProtassov et all |2002h and references therein for details. 



^In practice this was performed using the t clout simpars command in XSPEC 
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point made in section 13.1. 3[ that the Gaussian assumption is reasonable for these data. 



3.2.2. Automated fitting of GRB spectra 

Given the number of simulated datasets one must resort to an automated model fitting 
procedure. This has it self been the cause of some debate, with some authors (e.g §5 of 



Rutledge fc Sakd (|2003| )) claiming that automatic routines do not robustly find the best- 
fitting parameter values (minimum y^). The algorithm used by XSPEC for minimization 
is the Levenberg-Marquardt algorithm, which is efficient and very effective when the space 
is well-behaved (e.g. with only one local minimum). However, as this is a 'local' routine 
there is no guarantee of finding the 'global' minimum in x^, and it is possible that the results 
are biased by the presence of other local minima. For the present paper we have employed 
several additions to the standard Levenberg-Marquardt minimization algorithm in order to 
mitigate these problems. 

Once a local minimum in is found the surrounding region of parameter space is 
explored for signs of other minima. Each parameter in turn has its value increased and 
decreased until the is increased by at least Ax^ = 2.7, while simultaneously allowing 
the other parameters to vary in order to minimize Ax^. If any non-monotonicity in x^ is 
detected during this search the volume of parameter space explored is increased by increasing 
the value of Ax^. If during the course of this search Ax^ becomes negative (meaning there is 
a lower minimum nearby) the Levenberg-Marquardt algorithm is re-started from the position 
of this new minimum. The entire process is repeated by perturbing each parameter in this 
way until no further improvement can be made by the adjustment of any of them. 

The absorbed power law model (Mq) has only three parameters (photon index F, nor- 
malization and absorption column density), and in all cases finding the x^ minimum was 
straightforward using the above procedure. The alternative model Mi, which includes the 
emission line, required more care because the presence of a line with unknown energy may 
cause local x^ minima at different energies within the wide bandpass. An initial 'best guess' 
line energy was computed for each spectrum in the following way. An absorbed power law 
plus emission line model was constructed using the best-fitting parameters of model Mq and 
adding an unresolved emission line fixed at some trial energy Ei and varying the other pa- 
rameters (including the line normalization) to find the minimum x^- One hundred values 
of the trial energy Ei were used, evenly spread over the entire useful bandpass, and the 
value that recorded the lowest x^ was selected as the 'best guess' for the line energy. The 
enhanced Levenberg-Marquardt algorithm described above was then used to find the global 
X^ minimum starting from this position. 
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Simulation tests and comparison with interactive fitting demonstrated the automatic 
procedure described above was an efficient and very robust procedure for finding the global 
minimum. 



3.3. Rutledge and Sako Smoothing (RS) 



Rutledge &: Sakd (120031 ) proposed an alternative method for line detection using a 
'matched ffiter' to smooth the observed count spectrum with the aim of removing low sig- 
nificance noise and emphasizing any spectral features. The distribution of peak fluxes in the 
smoothed spectrum is then compared to the result of Monte Carlo simulations to calibrate 
their significance (p- value). 

The counts per PHA channel are extracted from the observed X-ray spectrum and 
then smoothed using an energy-dependent kernel ( a Gaussian having a FW HM equal to the 
spectral resolution of the detector; see equation 2 of iRutledge fc Sakd (120031 )) to produce the 
smoothed spectrum C{E). The distribution of C{E) is then calibrated using Monte Carlo 
simulations of spectra generated using the m ethod discussed in sect ion 13.21 that eniploys 
posterior predictive data sets. (We note that iRutledge fc Sakd l2003l and ISako et al. 



2005 



did not randomize the parameter values but used fixed MLE values to generate all their 
simulations. This is equivalent to assuming the posterior to be a delta function located at 
the best fit point, which is clearly a bad approximation in many cases.) Each simulation is in 
turn smoothed using the same energy kernel to produce C{E)sim,i- The C{E)sim values are 
then sorted in descending order for each PHA channel separately. Thus the 99*^ percentile 
limit of the C {E) sim,giobai is then found by extracting the 100*'' highest value of C{E)sim in 
each PHA channel. 

The smoothed observed spectrum, C{E), is then plotted alongside the n*^ percentile 
limits, which we have chosen for this analysis to be 90.00, 99.00, 99.90 and 99.99 %. Wherever 
C{E) exceeds a given limit then we have detected a 'feature' at that confidence limit. Thus 
a line would show up as a narrow excess whilst other thermal emission components will show 
up as broad excess, both of which are easily distinguishable. 



3.4. Comparison of the methods 

The three methods discussed above have different theoretical motivations, underlying 
assumptions and require different amounts of computing power. The Bayes factor method is 
based on a simple application of Bayes theorem combined with the Laplace approximation 
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and assumes uniform priors on the model parameters (or Jeffreys prior for the hne normal- 
ization). As discussed above, this may not be the optimal assignment. However, despite 
its possible drawbacks, the simple priors and Laplace approximation make the calculation 
extremely simple, requiring only the evaluation of equations [7] and [S] which require the values 
of and the covariance matrices for the best-fitting line and line-free models, and details 
of the free parameters and their allowed ranges. As such, it is useful as a 'quick and easy' 
test. The dependence of choice of priors may be assessed by comparing the results computed 
using the uniform and Jeffreys prior. 

By contrast, the RS and 'p'p'p methods require a large number of random datasets to be 
simulated and analyzed, and are therefore considerably more costly in terms of computing 
time. There is no compelling theoretical reason for applying a matched filter, as in the RS 
method, although it should be noted that the method, as implemented above, is calibrated 
using the appropriate posterior predictive distribution. The advantage of the RS method is 
that no model fitting is required, which is often a time-consuming process and can lead to 
biased results if not handled properly (section [3.2.21) . 



The vpy method is gro unded in the theory of Bayesian model checking (IGelman et al. 



19951 : iProtassov et al.ll2002l ) but requires time-consuming fits to be performed on each sim- 
ulated spectrum, and is therefore the most computationally demanding method by a clear 
margin. However, it is arguably the m ost rigorous in the sense that it is less sensi tive to the 
choice of priors than are Bayes factors (IGelman et al.lll996l : IProtassov et al.ll2002l ). and does 
not apply an ad hoc smoothing, as in the RS method, that may actually act to suppress real 
spectral features in some cases. 

The simulations used for both RS and 'p'p'p methods were generated assuming a Gaussian 
posterior for the three parameters of Mq, which, as discussed in section 13.1.31 is a good ap- 
proximation. Again, this approximation was made to increase computational efficiency, since 
Gaussian deviates are trivial to generate with the Cholesky method. In situations where the 
Gaussian approximation is not valid and/or the number of spectra is small enough that con- 
siderably more computing time may be spent on each, the vw method or Bayes factors may 



be co mputed using results from MCMC simulations (Ivan Dyk et al.l l200ll : IProtassov et al. 
2OO2I ) which allows for a more accurate evaluation of the posterior density. 



3.4- 1- Alternative approximate methods 



The statistics literature contains many methods developed for the purpose of model 
selection. In the introduction we listed four methods that have previously been applied 
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to the problem of line detection in X-ray data from GRBs. One method that has not, to 
our knowledge, been appl i ed sp ecificallv to GRB line detection is the Bayesian Information 
Criterion (BIG; ISchwarta (119781 ) ). This aims to approximate the logarithm of the integrated 
posterior probability for a model with k parameters given data with a sample size N. The 
BIG takes the form of the logarithm of the likelihood with a penalty term: 



BIC = - \n[p(D\e,M)] + {k/2)\nN 



(13) 



The model with the smallest BIG is favored. The difference between the BIG values 
for two competing models (often called the Schwartz criterion) is therefore S = — In A + 
(A/c/2) In (see footn ote EI) , and is a rough approximation to the logarithm of the Bayes 
factor (section 4.1.3 of iKass fc Rafterylll995l ). 



In the high count (large sample size) limit (see equation [5]) the Schwartz criterion be- 
comes 25* = — Ax^ + AfclnA^. Whether or not the BIG for model Mi is smaller than that 
for Mq is then equivalent to the criterion Ax^ > —AklnN. In the present case the data are 
selected with fixed N, and Ak = —2 for the addition of a fixed width line, such that the 
BIG is equivalent to applying the same Ax^ criterion to each spectrum, mechanically the 
same as the LRT, although with a different (generally higher) threshold value. Therefore, 
in the present context the application of the BIG would be equivalent to a slightly more 
conser vative application the LRT (see footnote fT5l) . However, as noted in iProtassov et al. 
(120021 ) and elsewhere, the BIG is often a poor approx imation to the integrated posterior 
probability, and as discussed by iKass &: Rafteryl (119951 ) is generally a worse approximation 
than the Laplace approximation employed to calculate Bayes Factors in section 13.1.21 



4. Results from an Iron line emitting source 



As a first demonstration of the above methods we applied them to a non-GRB Swift 
dataset. Ideally we would prefer to examine a source with a GRB-like spectrum, with a 
similar count rate, but containing a clearly identified emission line feature. However, it is 
difficult to find a source that meets all of these criteria. We chose the PG mode calibra- 
tion dataset (combining all av ailabl e data froni 10/05 /2005 to 02/09/2005) of PKS 0745-19 
JPe Grandi fc Molendj Jl999h and iGhen et Zl toO^ ). This test has some limitations as 
PKS 0745-19 is fainter than the GRBs analyzed in this paper and it is observed in a differ- 
ent mode. 



PKS 0745-19 is a galaxy cluster with a thermal spectrum and a known line at 6.07 keV 
in Swift's observations, which is a redshifted 6.7 keV iron line {z = 0.1). Even though the 
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underlying spectrum is thermal, with multiple temperature components, it can be modeled 
by an absorbed power law contii iuum where the po wer l aw index, F, i s 2.34]']qq3 [x^/i' = 



635/511). An additional mekal (IMewe et al.l (119851 ) and lArnaudI (119961 )) component, with 
kT = 0.191q'q5, slightly improved the fit with x'^/^ = 610/509. All spectral parameter errors 
are quoted at 90% confidence. 



Adding a Gaussian component to an absorbed power law fit naturally produced a signif- 
icantly improved fit to the data {x'^/u = 539/508) with T = 2.36t[{:[]^ and a line at Q.OYto^ol 
keV (width = 0.06]'^o;o3 keV). The spectral fit to this model can be seen in fig El This is 
supported by the Bayes factor of 7 x 10^^ for a single line being present. RS analysis of 
the spectrum, fig. [3], also clearly showed the presence of a Gaussian feature at ~ 6.07 keV 
with a significance far in excess of the 99.99% confidence limit. The ppp analysis placed a 
significance of > 99.99% on this feature. 

An interesting point to note is that there are shallow 'excesses' at ~ 0.6 keV and ~ 2.3 
keV in the RS plot (fig [3]), which are clearly not line features. Coherent, low level, positive 
excesses are also seen in the spectral fit at these energies (fig. [2]). Either the power law 
component is not modeling the data adequately at these points, the energy scale for this 
spectrum has an offset or the calibration files are less accurate around these two energies. 
Applying the gain f it function in XSPEC improves the model fits significantly by adding an 
offseto of -0.07 keV (no change to the slope). The absorbed power law model improves from 
X^ /v = 635/511 to 606/509 and the mekal component model improves from x^/^ = 610/509 
to 585/507. As a result the two shallow 'excesses' at ~ 0.6 keV and ~ 2.3 keV become far 
less prominent. 

The feature at ~ 0.6 keV could be attributed to the detector oxygen absorption edge 
at 0.54 keV. Applying the -0.07 keV offset brings the ~ 0.6 keV 'line' in conjunction with 
this edge, thus reducing its significance below the point at which we would consider it to be 
a real detection. We note that the ~ 2.3 keV feature is coincident in energy with the gold 
edge due to the XRT mirrors. We have confirmed that this feature is not due to any bad 
pixel or hot column issue 



^°http: / /swift. gsfc.nasa.gov/docs/heasarc/caldb/swift/docs/xrt/xrt_bias.pdf 
"http://swift.gsfc.nasa.gov/docs/heasarc/caldb/swift/docs/xrt/SWIFT-XRT-CALDB-01_v5.pdf 
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5. Testing the three methods and determining detection hmits. 

In this section we discuss the sensitivity limits of the three methods, i.e. the weakest 
lines that can be reliably detected with each of the three methods, for observations of the 
type discussed in section [2], of a 'typical' Swift era burst. This is done by simulating XRT 
data with a continuum spectral model typical of the GRBs observed with Swift, but including 
an emission line, and then applying the three methods described above for line detection. 

In order to generate the simulated data we use a fiducial spectral model comprising a 
power law with photon index F = 2.0, normalization (at 1 keV) of = 0.9 photons keV~^ 



s ^ an d an absorption column density of A^h = 1.8x10^^ cm ^ (see table 2 of ICampana et al. 



20064). These parameters are typical of the X-ray spectra of Swift era burst^lj. In order to 
measure the sensitivity of the three detection methods to lines in XRT data, spectral data 
were simulated using the above model plus one Gaussian emission line, and subjected to each 
of the three procedures. A range of values for line energy, normalization and intrinsic width 
were used in order to calibrate the dependence of the methods to the line parameter^. 

The ppp and RS method result in p-values with the conventional frequentist interpre- 
tation. If we set the detection threshold at a, and identify a detection as p < a then the 
rate of type I errors (i.e. false positive detections) will be a. For the purpose of sensi- 
tivity analysis we used a = 0.01, equivalent to a "99% significance" criterion. In contrast 
to these, the Bayes factor is the ratio of the marginal likelihoods of models Mi and Mq; 
in the case of uniform priors for the two models this is the ratio of posterior probabilities 
-^10 = p(^i|D)/p(Mo|D) where the probabilities are interpreted directly as probabilities for 
models Mq and Mi, respectively. 

For the purpose of numerical comparison with the p-values, the Bayes fact ors were 



conver ted into probabilities (assuming p(Mi|D) -|-p(Mo|D) = 1; see equation 3.19 of lGregory 



(120051 ) )■ and p(Mo|D) < 0.01 was taken as the criterion for detection. This is equivalent 



We have assumed a redshift z = for the fiducial burst spectrum. The average of the measured redshifts 



for Swift GRBs is higher than this (see http://www.astro.ku.dk/~pallja/GRBsample.html for the updated 



values). However, it should be noted that increasing z causes the effects of absorption by the host galaxy 
absorption (which tends to dominate the total absorption column) to shift out of the observed bandpass, 
meaning there is relatively more flux at lower energies (< 1 keV). The calculated detection limits should be 
representative of Swift era bursts although perhaps conservative at lower energies. 

The ranges of values used for the line simulations are as follows: Normalizations of 1 x 10^^ 100 
photons cm~^ s~^ taken in logarithmically increasing steps, line energies of 0.4, 0.6, 0.8, 1.0, 2.0, 3.0, 4.0, 
5.0, 7.0 and 9.0 keV, and intrinsic widths of 0.0 keV (i.e. unresolved), 0.2 keV (broad line) and 0.7 keV 
(broad continuum excess). 
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to p(Mi|D) > 0.99, and approximately equivalent to a Bayes fact or -Bm > 100, which i s 



conventionally taken as strong evidence in favor of Mi over Mq (iKass fc Rafteryl Il995l ). 
However, we stress that the interpretation of p-values and Bayes factors are fundamentally 
different. A p-value is the tail area of the probability density function of the test statistic, 
assuming a null hypothesis (Mq) is true, and is used to decide whether or not to reject 
the hypothesis. As such, a p- value is not the probability for the model Mq, instead it 
corresponds to the frequency of more extreme test statistics (e.g. Ax^) given a large number 
of repeat experiments (assuming the null hypothesis). By contrast, p(Mo|D) is the posterior 
probability for model Mq based on data D and the priors (in the present case we used an 
approximation thereof), as p(Mi|D) is for Mi, and Bayes factors are used to select between 
two models based on the ratio of these two. This fundemental difference in the interpretation 
of Bayes factors means there is no expectation that a is the frequency of type I errors from 
a large number of repeat observations when using a p(Mo|D) < a criterion. 

In § 13.1.31 we confirmed that using the Laplace approximation assumption in the calcu- 
lation of the Bayes factor was valid for the fiducial absorbed power law spectral model. The 
same was also found to be true of the spectra with simulated Gaussian lines at, and above, 
the detection limit detailed above. 

For each value of the line normalization we calculated the Bayes factors for 50 inde- 
pendent simulations and calculated the p(Mo|D) values for each. We then averaged the 
p(Mo|D) values at each normalization and linearly interpolated between points at adjacent 
normalization values to map Pq as a function of normalization. The limiting sensitivity was 
taken to be the normalization at which the mean p(Mo|D) value falls below 0.01. 

Figure m shows the detection limits for an intrinsically narrow line {W = 0) at different 
energies for spectra with ~ 800 and ~ 1600 counts (left and right panels, respectively). 
The limiting sensitivities are shown in units of equivalent width (keV), which is easier to 
interpret physically, than the absolute normalization, by comparing the normalization to the 
underlying continuum model. Figures [5] and [6] show the detection limits for different line 
widths [W = 0.2 and 0.7 keV, respectively). The Bayes factor points in these figures have 
been calculated using the uniform prior, rather than the Jeffreys prior. See § l6.11l for further 
discussion on the effect of using the two different priors in the calcualtion of the Bayes factors 
for the observed data sets. 

The ppp and RS methods require a large number of spectral simulations in order to 
calibrate their distribution and estimate the p-value. The computational demands of thir^ 



^''To give a specific example, for the simulation and fitting methods described in sections 13.2.11 and 13.2.21 
a set of iV = 10^ simulations takes ^ 1 day on a top-range PC. 
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are such that it was not feasible to produce a sufficiently large set of simulations to carry out 
the methods on each and every line spectrum (e.g. which includes several spectra at each trial 
value of line energy, width and normalization, for both ~ 800 and ~ 1600 count spectra). 
We therefore constructed two libraries of 10^ simulations, one for ~ 800 and one for ~ 1600 
count spectra, that could be used for each test. These were constructed by simulating an 
appropriate spectrum based on the fiducial model, and using this to generate the posterior 
predictive distribution from which to draw 10^ simulations following the recipe discussed in 
section [3.1.21 These libraries were then used to calibrate the distribution of the Ax^ statistic 
for the ppp method and thus to calculate the value of Ax^ that corresponds to a p-value of 
0.01. Similarly these libraries were used to compute the 99.00% significance contour from 
the fiducial model for the RS method. We point out here that these simulation libraries were 
used only for the purposes of comparing the different algorithms. For the analysis of real 
observations (discussed below), each observation was assessed using independently generated 
simulations matching the particular observational parameters. 

For the ppp method each of the spectra containing a line was fitted with an absorbed 
power law with and without an additional Gaussian component, and the change in noted. 
The Ax^ values were averaged at each normalization, and these points linearly interpolated, 
to map the Ax^ as a function of normalization. As with the Bayes factor, the limiting 
sensitivity was taken to be the normalization at which the mean p-value falls below 0.01, 
calculated using the appropriated value of Ax^ value from each simulation library. The 
limiting sensitivity as a function of energy is shown as green dotted curves in Figures IH O 
and [HI for different configurations of line parameters. 

For the RS method each line spectrum was smoothed individually. The C{E)sim values 
over an energy channel range equal to the central energy, Eune, ± line width were extracted. 
These values were compared to the 99.0% confidence limit over the same energy channel 
range found from the appropriate simulation library. The number of channels within this 
range where C{E)sim > C{E)gQX) was recorded for each simulation. The detection limit was 
taken to be the lowest line normalization where N{C{E)sim > C{E)gg,o) = 0. 

Analysis of the (line free) library simulations showed the Bayes factors produced < 1% 
false positives when Biq > 100 was used as a detection criterion. This shows the method 
is, if anything, slightly conservative as expected given the conservative Amax assumption. 
Conversely, the false negative detection rate is negligible above the detection limit. 

As expected the detection limits are higher for the spectrum with a lower number of 
counts, by a factor of ~1.5. For the fiducial spectral model used here the optimum energy 
range for detecting lines is 0.4 — 6 keV, where the line only requires a contribution of a few 
% of the total spectral flux. In the best cases (1600 counts and narrow line) a line with an 
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equivalent with as small as ~ 40 eV may be detected around 1 keV (observed frame) at 
99% significance (in a single trial), whereas only very strong lines may be detected between 
6 — 10 keV. Additional simulations were carrie d out with a highe r absor ption column density 
(1.0 X 10^^ cm~^; the mean values stated in iReichart &: Pried (|2002| ) assuming that long 
bursts occur in molecular clouds). The dependence of line detection with respect to energy 
for all three methods were the same at energies >1 keV. However we note that simulating 
the spectra with much larger absorption columns significantly degraded the ability to detect 
lines features below 1 keV. 



6. Results from Swift archival GRB afterglow data 

Our sample covers a subset of 40 GRBs, out of the total of 153 from GRB 050128 up to 
GRB 060510B, which were selected for the quality of their WT mode data (see §[2]). Some 
bursts only contained sufficient data for a single WT mode spectrum to be analyzed, whilst 
the majority contained sufficient data to be time-sliced into multiple spectra (see § [2])- In 
total 332 spectra were analyzed. We sample a range of energies and time spans even though 
the complete redshift distribution of this data set is unknown. The subset of this sample 
with known z indicates that we are typically probing the region between T+0 s to T+~500 
s (or up to T+few ks if the burst is very bright) post trigger and between ~ 1.0 and ~ 50 
keV in the rest frame of the burst. Throughout this section error bars indicate nominal 90% 
confidence limits on one interesting parameter. 

All the data were fitted using automated procedure described in section 13.2.21 and the 
solutions checked by hand. In practice four models were fitted to each spectrum: (1) absorbed 
power law; (2) absorbed power law plus unresolved Gaussian emission line; (3) absorbed 
power law plus variable-width fine; (4) absorbed power law plus blackbody. The results 
presented below focus on the line models, and we found that the blackbody parameters were 
in general very poorly constrained. Ideally, we would like to apply all three methods to all 
332 WT spectra to assess the significance of lines (or other) features in the data. But, as 
discussed above, the RS and especially ppp methods are computationally demanding and so 
it was not practical to apply these methods to every spectrum. 

The (approximate) Bayes factor method, being computationally economical, was applied 
to every spectrum, while the more computationally expensive RS and ppp methods were 
applied only to subsets of the data. In particular, any spectrum that showed a Bayes factor 
Bio > 1 (in favor of a line), or a Ax^ > 4.61 upon inclusion of a linfl was considered 



i5/\t^2 _ 4.60517 is the 90th percentile for the distribution with 2 degrees of freedom (jPress et al. 
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for more detailed analysis. These were deliberately chosen to be extremely relaxed selection 
criteria (especially so given the large number of independent tests, see below), so as to avoid 
removing any plausible line candidates and only remove those spectra without any hint of a 
line, and to counteract the conservative nature of the Bayes factors (section [S]). Indeed, this 
screening effectively reduced by a factor ~ 4 the sample of spectra worth considering in more 
detail. We re-iterate that no judgement about the presence/absence of a line in a spectrum 
was made purely on the basis of the Bayes factor method, which, as discussed above, is an 
approximation and is sensitive to the choice of priors. Only spectra with a low Bayes factor 
(-Bio < 1) one? little improvement in the fit statistic upon including a line (Ax^ < 4.61) were 
not considered for further analysis. This subsample was then subjected to the RS method 
with a low detection threshold {p < 0.1, i.e. a 90% single trial significance, again very weak 
given the multiple trial effect). This further reduced the sample size to a level where the 
rigorous but computationally expensive ppp method could be applied. 

As stated above, this screening was only necessary to reduce the sample to a manageable 
size for ppp analysis. Numerical tests showed that the ppp method invariably gave a higher 
p-value (i.e. lower significance) than the RS method, and so no data that might have shown 
a detection with the rigorous ppp method would have been lost by the selection process. 

The large number of spectra examined means the effects of multiple trial must be in- 
cluded in the analysis. For example, to reach a global detection significance of only 90.0% 
given a sample of 332 spectra, we would require a single trial significancJ^ in excess of 
99.97%. Of the 332 spectra from 40 GRBs, 12 spectra from 10 GRBs gave a single trial 
detection of > 99.9% significance in at least one of the methods. As the best line candidates 
in the sample, we now consider each of these in turn. (All significances are single trial values, 
unless otherwise stated.) 



19921 ). As such, it corresponds to a 90% detection 'significance' {p < 0.10) in a classical likelihood ratio test 



(LRT) when including two additional parameters (see footnote [8]) . T he LRT should not be used directly for 
the purposes of detecting an emission Hue (for reasons discussed in lProtassov et al1l2002h . but in practice 
the p- value calculated from the analytical test is usually within an order of magnitude of the value calibrated 
using the ppp method. It is therefore extremely unlikely that a dataset producing Ax^ < 4.61 would yield 
a solid detection (e.g. p < 10"^) after ppp analysis. 

^^Calculated using the standard Bonfcrroni-type correction factor: pi = 1 — (1 — pnY^'^ , where pi is the 
single trial p- value that gives pN as the rate of type I errors in a set of N independent trials. This sometimes 
known as the Sidak equation. In this limit of small pn and large N this tends to pi = p^/N. 
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6.1. GRB 050730 



A single Gaussian feature was detected in the spectrum e xtracted from 



T+792s, which was concurrent with a flare in the WT mode data (IStarling et al.l (120051 ) and 



'+692s to 



Pandey et al.l (120061 ) ). An absorbed power law plus a broad Gaussian {a = 0.34_Q;^g keV) at 



1.14to;^4 keV provided the best fit to the data with xV'^ = 47/52 (table [I]). When the line 



width was restricted to below the detector resolution a Gaussian feature at 0.73]';o;o3 keV 
was detected {x^ = 57/53). 

The Bayes factor was Biq = 300, favoring a line. The RS method (fig. [7]) indicated that 
a line is present in the spectrum at ~ 0.7 keV with a confidence of 99.90%. This compares 
favourably to the parameters found in the spectral fit when the Gaussian width was restricted 
to a value below the instrumental resolution. There is no evidence for the broader feature 
found when the width of the Gaussian was a free parameter (see inset to fig. [7]). 

A ppp analysis was carried out in both cases. The significance of the unresolved-width 
and free-width Gaussian features were found to be 88.50% and 99.92% respectively. It is 
surprising that the ppp analysis appears to favour the wider line at i^une = l-14^o;44 keV, as 
there is no evidence of a feature with this energy in the RS plot. However, we note that the 
large errors on this line energy are consistent with a feature at 0.73loo3 keV at the limit of 
their range. 

Applying the gain fit function to this spectrum resulted in an improved fit [Ax^ = 6) 
for an unresolved-width line feature at 0.68]!^g o4 keV, with an offset of -55 eV (all other 
spectral parameters were unchanged within previous limits). Combining this energy offset 
with the error on the line energy is not sufficient to prove an association with the oxygen 
absorption edge. Applying the gain fit function to the free- width Gaussian model was 
inconclusive, with regards to an association to the oxygen absorption feature, owing to the 
poorly constrained line energy of < 1.00 keV. (For further discussion on the application of 
the gain fit function to this and other GRBs see §7.) 



Th e redshift for this burst was re ported as z = 3. 967 (Chen et al.l (120051 ) . iHolman et al. 



(120051 ) ■ iProchaska et al.l (120051 ) and IStarling et al.l (120051 )). Further fits were conducted 
wi th two Njj columns origin ating from the Galactic column (wabs, fixed at the value given 
by lDickey &: LockmanI (Il990l )) and the host galaxy (zwabs). This had the effect of marginally 
improving the fit for the absorbed power law model (x^/ = 61/55, Ax^ = 3) with a Galactic 
column density fixed at 3.21 x 10^° cm~^ and a host galaxy component of 9.80lg2o ^ ^^^^ 
cm~^. All of the other spectral parameters were the same as the previous fit within the 
limits. The fit to the other models, containing Gaussian components, did not change sig- 
nificantly and the parameter values were the same within the error limits. Bayes factor 
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analysis including the zwabs component indicated marginal evidence for line being present 
(-Bio = 5). Applying the additional zwabs component to the RS method (see fig [8]) also 
decreased the significance of the 0.73 keV feature from 99.90% confidence (dotted line) to 
99.0% (solid line). A ppp analysis, taking the zwabs component into account, found that the 
significance of the free-width feature had decreased to 99.49% (i.e. ~ 2.8a detection) in this 
single trial. We conclude that the line detection (unresolved or free-width) in GRB 050730 
is not significant at 3a, and note that the redshift-corrected line energy does not correspond 
to a K-shell transition of a common element. 



6.2. GRB 060109 

This burst had insufficient fiux to produce multiple spectra therefore we considered 
the dataset as a whole. The spectrum covers data from T-|-109s to T+199s. An absorbed 
power law plus a narrow Gaussian at 0.741q;q3 keV (width restricted to below the detector 
resolution, x^/^ = 40/40) and a free-width Gaussian at < 0.72 keV (width = 0.23+o;og keV, 
X^/^ = 35/39) were equally good fits to the data (tabled]). 

The Bayes factor for the unresolved-width Gaussian model indicated the presence of 
a line {Bio = 200), however the same analysis on the free- width Gaussian was much less 
convincing {Biq = 3). The RS method indicated that there may be a feature at ~ 0.7 keV 
with a significance of 99.90% (see fig [9]). However, the ppp method gave only 88.99% and 
99.28% significance for unresolved (fixed) and free-width Gaussian lines, respectively. In §6.1, 
we showed that the significance of a similar feature decreased below 3a when the spectral fit 
was changed to include an absorption component at the redshift of the host galaxy. We will 
show that this is generally true for those GRBs for which a redshift is known. Unfortunately, 
in this case, the redshift is not known and we cannot determine whether or not the same is 
true. 



6.3. GRB 060111A 

The data from this burst were split into 11 spectra, covering several fiaring events that 
showed significant spectral variation during the observation. The Bayes factor {Bio = 0.05) 
gave no evidence for a free-width line (at 0.65lo;og keV, W < 0.13 keV) in the spectrum 
covering T+174s to T+234s, despite it producing a modest improvement in the fit (Ax^ = 8; 
table [I]). The RS results {figM), a feature at ~ 0.65 keV with > 99.90% confidence. A 
further feature at 0.79l;g;Q^ keV {W < 0.15 keV) was detected in the spectrum covering 
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T+319s to T+339s. The Bayes factor indicated that the presence of a hne in the second 
spectrum was unhkely {Biq = 0.1) but the RS method (fig [11]) suggested an additional 
spectral feature. 

Whilst the ~ 0.65 keV feature for T+174s to T+234s and the ~ 0.79 keV feature in the 
T+319s to T+339s both look promising from the RS method, the ppp analysis showed that 
they were only 85.13% and 99.56% significant, respectively, not strong detections given the 
number of trials (see above). 



6.4. GRB 060115 



This burst had insufficient fiux to produce multiple spectra; therefore we considered 
the dataset as a whole. The spectrum covers data from T+121s to T+253s. An absorbed 
power law plus a Gaussian at 0.811q;o7 keV with a width of O.IOIq^os keV provided the best 
fit to the data with x^/^ = 82/77 (tabled]). The Bayes factor gave no evidence for a line 
{BiQ = 0.03), but the RS results (fig. [12]) indicated that there was a feature at ~ 0.75 keV 
at the 99.90% significance. However, ppp analysis gave only 96.16% significance. 



A redshift of z = 3.53 was reported by jPiranomonte et al.l (120061 ). Further fits were 
conducted with two ./Vh columns or i ginat ing from the Galactic column (wabs, fixed at the 
value given by iDickey fc LockmanI (jl990l )) and the host galaxy (zwabs). This led to no 
change in the statistical fit nor parameter values for an absorbed power law model or models 
containing Gaussian components. We conclude that the line detection in GRB 060115 is 
only moderately significance in a single trial, and not significant (to 30") in multiple trials, 
and note that the redshift-corrected line energy does not correspond to a K-shell transition 
of a common element. 



6.5. GRB 060124 

A precursor ~ 570 s before the main burst peak allowed Swift's nar row-field instrument s 



to be positioned on the GRB location ~ 350 s before the burst occurred ([Romano et al.ll2006l ). 
Therefore the WT mode data covered both the prompt emission from the burst as well as 
a portion of the afterglow phase. The fiux detected over the observation was sufficient to 
produce a time series containing 46 spectra. Bayes factor analysis indicated that eight of 
these showed evidence for additional spectral features and a further nine showed evidence 
from the raw Ax^ improvements. However, RS and ppp analyses carried out on all of these 
potential line spectra revealed only one with an acceptable detection (with both methods 
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giving a significance of > 99.90%). Tliis spectrum spanned T+537s to T+542s (i.e. occurring 
just prior to tlie main burst peak). Tfie best fit model to tliis spectrum was an absorbed 
power law plus a broad (a = 0.48toii keV) Gaussian component at 2.30l^o23 keV (tabled]). 

The Bayes factor was Biq = 20 for a free-width Gaussian feature in this spectrum. RS 
results (fig. [13]) showed a 99.99% significance feature at ~ 2.55 keV. A ppp analysis indicates 
that the feature is significant to 99.97%. 

Whilst this appears to be a significant detection it seems to be very broad for a sin- 
gle li ne feature, requiring a v elocity dispe rsion of the order 0.5c. Using the redshifts of 
0.82 jMirabal fc Halpernlliooel ) and 2.297 JCenko et al.ll2006h it is possible to identify this 



feature with Ka emission of Calcium (4.10 keV) or Cobalt (7.5 keV) respectively. It could in 
principle be a series of unresolved line features, a thermal component or indicating a break 
in the spectrum. Fitting the spectrum with a blackbody component (kT = 0.76^°;]^^ keV) 
did not provide a good fit {x^/v = 61/45) nor does an absorbed broken power law model 
(xV^ = 56/45). 

A further possibility is that it could be due to a poor fit to the gold M-edge as seen 
in §11 However, applying an energy offset to the data does not significantly improve the 
absorbed power law model (Ax^/z/ = 73/45, offset = -0.08 keV, no change to the slope). 



6.6. GRB 060202 

This burst contained sufficient flux to extract 18 spectra. Of these spectra only one, 
spanning T+429 s to T+529 s, appeared to contain an additional spectral feature. The 
Bayes factor was Bio = 300 in favor of a single free-width line. An absorbed power law 
plus a broad (< 0.34 keV) Gaussian feature at 0.94+°;°^ keV {x^/p = 96/100) was a slightly 
better fit than an absorbed power law alone (x^/ = 109/103; see tabled]). RS results for the 
T-l-429 s to T+529 s data (figdl]) indicated a broad feature at ~ 0.95 keV, which exceeds the 
99.99% confidence interval. A ppp analysis of the same data places a significance of 99.74% 
on this broad feature. No redshift value has been reported for this burst thus we were unable 
to perform a well constrained two component absorption fit. 



6.7. GRB 060210 

This burst contained sufficient flux to extract a time series containing 8 spectra. Of these 
spectra only one, spanning T+233s to T+353s, appeared to contain an additional spectral 
feature. A model containing a Gaussian feature at 0.671o;q4 keV (width = 0.06^g;o3 keV) 
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was a much better fit than an absorbed power law alone (Ax^ = 13; table [T]). The Bayes 
factor was Biq = 1. The RS results (fig [T5l) indicated that a feature at ~ 0.65 keV with a 
significance of 99.99%. A ppp analysis of the spectrum indicated that the same feature is 
significant to 99.83%. 



A redshift of 3.91 was reported by ICucchiara et al.l (120061 ) for this burst. A two com 



ponent absorption fit was carried out on the data. This produced a significantly improved 
fit to the absorbed power law (T = 2.5011q \1) model with a host A^h column contribution of 
5.65l|^j^ X 10^2 and xV^^ = §0/72 (Ax^ = 18 compared to the fit with free Galactic ab- 
sorption only). Adding a Gaussian component to this gave T = 2A6~^q\1, a host absorption 
column of S.Tll^Q gg x 10^^ cm~^ and x^/^ = 76/69 (Ax^ = 5 compared to the fit with free 
Galactic absorption only). The addition of the zwabs component did not change the energy 
of the feature but was only able to place an upper limit of < 0.10 keV on its width. Bayes 
factor analysis after allowing for a zwabs component indicated no evidence for an additional 
spectral feature {Biq = 5 x 10""^). We can conclude that this feature is most likely a false 
positive detection. 



6.8. GRB 060218 



Campana et al.l (l2006bl ) have reported on the association of this burst with SN2006aj 



and the presence of a thermal component in the X-ray spectrum in great detail. Our analysis 
concurs with their results. The data were split into 53 time intervals, from which the Bayes 
factor analysis indicated additional component in the spectrum in all data from ~T-|-750 s 
(with BiQ > 50). This was confirmed by RS and ppp analysis. The RS results (fig[IS], T-l-159 
s to T+2770 s) indicated that this feature is unlikely to be a Gaussian emission line as its 
profile was too broad. It is possible that it could be a series of unresolved lines, however, 
a power law plus blackbody component gave the best fit to all of the spectra suggesting an 
additional spectral feature. Similarly, individual time-slices (see figHZl T-|- 2359 s to T+2409 
s, for one such example) show th e presence of this broad feature, which appears to evolve 
over time (jCampana et al.ll2006bl ). 



6.9. GRB 060418 

A time series of 12 spectra were extracted from this GRB, two of which appear to contain 
additional spectral components. These were the spectra spanning T+119 s to T+129 s and 
T+169 s to T+194 s. 
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A Gaussian component at 2. 42^*^003 keV improved the fit to the T+119 s to T+169 s 
data by Ax^ = 16 (see tabled]), although the Bayes factor was unconvincing {Biq = 0.05). 
The RS analysis (fig. [TS]) showed a feature at this energy that clearly exceeded the 99.99% 
confidence limit. A ppp analysis found 99.85% significance for the same feature. However, 
as noted previously in the analysis for GRB 060124 and PKS0745-19 (§!]), a feature at this 
energy is coincident with the gold M-edge. 

A similar improvement in the fit was found for the second spectrum (T+169 s to T+194 
s), with Ax^ = 10, the Bayes factor was more promising {Biq = 30). An unresolved- width 
Gaussian at < 0.75 keV provided the best fit to this spectrum with x^/^ = 43/49 (tabled]). 
RS (fig. dH]) and ppp analysis supported the presence of this feature at the 99.99% and 99.98% 
confidence limit respectively. 

The 2.42 keV feature of the T+119 s to T+129 s spectrum can be explained by the gold 
M-edge but the 0.69 keV feature of the T+169 s to T+194 s spectrum cannot be matched to 
another elemental absorption edge in the same manner. Two component absorption fits were 
carried out with a A^h column densit y of 9. 17 x 10^ ° cm~^ from our Galaxy and a contributio n 



from the host galaxy at ^ = 1.49 (IDupree et al.l (120061 ) and IVreeswijk fc JaunsenI (120061 )). 



This produced a significant improvement in the absorbed power law model fit, which gave 
r = 2A8toil and a host A'h = 0.73l[J:| x 10^2 cm'^ {x^/u = 53/52, Ax^ = 9 compared 
to the fit with free Galactic absorption only). The addition of a zwabs component to the 
Gaussian model gave T = 2.22lgj5, a line with an energy of < 0.65 keV and width of 0.47]'ro;og 
keV and a host absorption < 0.21 x 10^^ cm"^ {x^/^ = 46/49). The Bayes factor for the 
spectra containing the zwabs component indicates that the odds of an additional spectral 
component have been significantly reduced to -Bio = 1- We can conclude that this feature 
is most likely not real, but a spurious detection due to the baseline assumption of no host 
galaxy absorption. 



6.10. GRB 060428B 

Data from this burst were split into two sets, T+212 s to T+252 s and T+252 s to T+418 
s. An absorbed power law model was a poor fit to the first spectrum with x^ /u = 78/63 whilst 
an absorbed power law plus a Gaussian feature at 0.76lg;og keV (width = 0.09]'^o;q3 keV) was 
a much better fit with x"^/^ = 63/60. However, the Bayes factor was less encouraging with 
BiQ = 0.1 against a line feature. The RS analysis (fig. [20]) indicated the presence of two 
possible features; one at ~ 0.75 keV at a significance of 99.99% and another at ~ 0.90 keV 
at 99.90%. However, no stable spectral fit could be found using an emission line at ~ 0.9 
keV, hence it was not possible to calculate a Bayes factor, nor calculate the Ax^ needed for 
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a ppp calculation. A ppp analysis of the feature at 0.7Qt.om yielded a significance of 
99.85%. 

The second spectrum, T+252 s to T+418 s, was best fitted by an absorbed power 
law model {x^/ty = 58/64, tabled]) and the Bayes factor gave only very weak evidence to 
indicate a line (-Bio = 3). RS analysis (fig [21]) indicated a possible feature at ~ 0.7 keV with 
a significance of 99.90%, however a ppp analysis placed a much lower significance of 95.87% 
on this. 

No redshift value has been reported for this burst, preventing us from performing a 
constrained two component absorption fit. This could potentially determine if the features 
at ~ 0.7 keV are due to poor modeling of the absorption continuum due to not including a 
component from the host galaxy. 

6.11. Use of alternative prior 

In section [3.1.41 we discussed two different choices for assigning an uninformative prior to 
the line normalization. The approximate Bayes factors given above were calculated assuming 
a uniform prior for the line normalization, but using the Jeffreys prior did not change the 
results significantly. For example, GRB 060115 changed from Biq = 0.03 with the uniform 
prior to Bio = 0.05 with the Jeffrey's prior. At the other extreme, the favorable Bayes 
factor of Bio = 300 for GRB 060202 (T+ [429-529] s) using the uniform prior was virtually 
unchanged. Spectra that were not included for further analysis, due to a low Bayes factor 
{Bio < 1-0) S'lid a small improvement (A^^ < 4.61) were similarly affected by a change in 
priors. In general the Bayes factors changed very little between uniform and Jeffreys priors, 
reflecting the fact that typical best-fitting line normalizations were usually ~ 10% of the 
total flux (see section [3.1.4p . 

7. Discussion of Swift XRT results 

The previous section shows that of 332 WT mode spectra analyzed by our methods only 
12 produced possible detections at > 99.90% (single trial). These detections were tightly 
clustered around two energies in the observer frame: 0.64-0.94 kev (10/12, figs [22] and [23]) and 
2.30-2.49 keV (2/12, fig [21]), with equivalent widths of ~0.9 keV and ~0.5 keV respectively. 

The coincidence of many spectral feature detections close to 0.7 keV is suspicious as 
we would expect intrinsic GRB emission line features to be located at different observed 
energies, as the GRBs span a large range of redshifts. This clustering strongly hints at 
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an instrumental origin. Modeling the WT mode spectra with an energy offset, in case of 
imperfect bias subtraction at the processing stage, improved the fit statistics for an absorbed 
power law model (average Ax^ ~ 5). However, even if the combined error on the line energy 
and the offset corrections (1 to 70 eV) are taken, this would still not be enough to provide 
a plausible association with the oxygen K-edge as seen in the PKS 0745-19 example (§4). 

An additional absorption component, at the host galaxy redshift, was applied to those 
candidates with a known redshift measurement. In every case the feature at 0.7 keV became 
insignificant and we expect that the same reduction in significance would occur if we were 
able to conduct well constrained two absorption component fits to the GRBs with unknown 
redshifts. It could be argued that this decrease in line significance stems only from the 
increased complexity in the model. If this were the case we should see an overall increase in 
the level of the RS contours over the whole energy range. In fig|H]we can see that the effect 
of the adding a zwabs component is not uniform across the energy range. It has negligible 
effect at energies >1.2 keV. Below this energy the zwabs component acts to increase the 
total absorption at very low energies (<0.55 keV) and decrease it in the 0.55 keV to 1.2 keV 
range. Thus the addition of the second absorption component is imposing a real, energy 
dependant effect on the confidence contours, rather than increasing them uniformly across 
the whole energy range. We conclude that the absorption is not being modeled accurately 
at low energies by assuming that all of the A^h column is at a redshift of zero. 

We have confirmed that the features at 2.3 keV are not due to bad pixel or hot col- 
umn issues. The two detections have the following single-trial significances: 99.97% (GRB 
060124) and 99.85% (GRB 060418). Taken in the context of all the trials performed these 
significances become 90.09% (< 1.7cr) and 50.20% (< 0.7cr) respectively. There is no signifi- 
cant impovement to the model fits if the gain fit function is applied. Adding a blackbody 
component to the underlying absorbed power law or using an absorbed broken power law 
does not yield a significantly improved fit either. We note, however, that both of these 
features are coincident with the gold M-edge complex. 

Since all of the features are found to be narrowly clustered around two energies in the 
observer frame, one of which was also found in the PKS 0745-19 example, it is our conclusion 
that none are real detections of emission lines in GRB spectra, but are instead either due 
to residual calibration issues, imperfect bias subtraction at the processing stage or incorrect 
modeling of the host galaxy absorption column (most likely for the 0.7 keV features). 

In addition we would expect of the order of 3 false positive detections at 99.00% for 332 
spectral analyses. 

It is also interesting that the majority of GRBs with potential features occur in the first 



- 33 - 



few months of 2006. There are no physical, instrumental or calibration issues associated 
with that period of operations that could explain such temporal clustering. However, we 
note that the actual response of the XRT CCD is possibly evolving while the response is 
modeled by the calibration files as being constant. 



Butlerl ( 120071 ) has recently published analysis citing the detection o f line complexes in 



GRBs 050714B, 050822, 060202 and 060218. However, the same paper and lButler fc Kocevski 



(120071 ) also provided alternate reasons for these apparent lines detections (thermal compo- 
nents or broken power laws). ORB 050714B was not included in our selection of bursts as 
the WT mode spectrum did not contain sufficient counts to meet our minimum criteria. 

We found no compelling evidence from GRB 050822 using any of our analysis methods 
to suggest th at there were any line features in these data, hence it does not feature in ^ 



Butlerl ( 120071 ) quotes a significance of AAa for a complex of 5 lines (0.81, 0.91, 1.04, 1.23 
and 1.49 keV) for one spectrum spanning T+489.5 s to T+509.4 s, with equivalent widths 
of 82, 142, 194, 221 and 265 eV respectively. We note that the features are spaced 100 eV 
apart, though no obvious physical explanation for this presents itself. It should be noted that 
Butler's spectra (~ 500 counts) is contained within our analysis from T+471 s to T+661 
s (~ 860 counts); our spectrum cannot be sub-divided further whilst still being directly 
comparable to the rest of the data analysis presented in this paper. None of the power law 
(or cutoff power law) models tested provide a good fit to the data (table [2]), though cutoff 
power laws appear to give a much better fit to the data. Adding a blackbody component to 
an absorbed power law gave poorly constrained parameters: T < 1.98; kT = 0.171q;q2; and 
Nh = 4.58tlfs X lO^o cm-2 {x^u = 46/38). 

RS analysis of our spectra compared to a base model of an absorbed cutoff power law 
shows a feature at ~ 0.7 keV at > 99.99% significance. However, ppp analysis only placed 
a significance of 99.74% on the same feature. In addition the appearance of a feature at 0.7 
keV indicates that the absorption column may not have been modeled accurately, as seen 
in other bursts in the previous section. As there is no report of a redshift for this burst we 
cannot confirm this by carrying out a well constrained two component absorption fit. 

060202 and 060218 have both been discusse d in the previous sectio n. We agree with the 



presence of a blackbody component in 060218 (ICampana et al.ll2006bl ). The single feature 
found in 060202 only occurred in 1 of the 18 time sliced spectra with a significance of 99.74% 
(single trial). If we consider the multiple trials carried out this drops to a significance of 
95.4%. Our analysis indicates that it is a broad feature {a < 0.34 keV) at 0.9 keV rather 
than a series of reso lved or overlapping lines, however it does occur at the same energy over 



which iButlerl (120071 ) reports 4 individual narrow lines. We find no evidence for the reported 



feature at 4.70 ± 0.07 keV. 
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8. Conclusions 

Analysis of the galaxy cluster PKS 0745-19, which has a known 6.07 keV emission line, 
produced a convincing detection by all methods and also uncovered two further features at 
0.6 keV and 2.3 keV. Both features are hkely to be due to an energy scale offset that causes 
the instrumental oxygen and gold absorption edges respectively to be poorly fitted (§|1]). A 
series of simulations over a range of emission line parameters has allowed us to estimate the 
sensitivity to Gaussian-like features, both broad and narrow. For all three methods, using 
GRB parameters typical for Swift bursts, the optimum range for emission line detection was 
found to be 0.4 — 6 keV, with line equivalent widths as low as ~ 50 eV detectable in principle 
from data with only ~ 1600 counts. 

Of the 332 WT mode spectra from real GRBs, only 12 produced possible detections at 
>99.90% (single trial). These were all located around two energies in the observer frame: 
0.7 keV (10/12) and 2.3 keV (2/12). The coincidence of many spectral features close to 0.7 
keV is suspicious as we would expect intrinsic GRB emission line features to be located at 
different observed energies, as the GRBs span a large range of redshifts. For those candidates 
with a redshift measurement the feature at 0.7 keV becomes insignificant once an absorption 
component at the redshift of the host galaxy is apphed. We expect that the same reduction in 
significance would occur if we were able to apply well constrained two absorption component 
fits to the GRBs with unknown redshifts. The 2.3 keV features are thought to be associated 
with the gold M-edge. 

Since all of the features are found narrowly clustered around two energies, one of which 
was also found in the PKS 0745-19 spectrum (2.3 keV), it is our conclusion that all of these 
features are either due to calibration issues, imperfect bias subtraction at the processing 
stage or incorrect modeling of the host absorption column (most likely case for the 0.7 
keV features), rather than GRB emission line detections. The only non-power law emission 
component we accept as intrinsic is the blackbody component detected in GRB 060218. 
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gratefully acknowledges support from a PPARC studentship. SV, JPO, ER, KLP, AN, DC, 
MRC, PE and RS acknowledge PPARC support. This research has made use of data obtained 
through the UK Swift Data Archive, provided by the University of Leicester. 




Fig. 1. — Marginal posterior distributions for the continuum parameters of an absorbed 
power law fit to an XRT spectrum of GRB 060124. The contours enclose 80, 50, 20, 10 and 
5% of the distribution (and therefore correspond to 20, 50, 80, 90, 95% credible regions for the 
parameters). The filled contours were computed assuming the posterior is a Gaussian. The 



hollow contours were computed using MCMC simulations from the routine of Ivan Dyk et al. 



(I2OOII ). The two distributions are clearly very similar. See section [3.1.31 for details. 
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Energy (keV) 



5 



Fig. 2. — Spectral fit to PKS 0745-19 with an absorbed power law plus a narrow Gaussian 
emission line model. The redshifted Iron line at 6.07 keV is clearly visible. Note also the 
residuals at 0.6 keV and 2.3 keV, which are thought to be due to the oxygen and gold edges 
respectively. 
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Fig. 3. — PKS 0745-19: RS method. Confidence contours mark the significance of the 
spectral features. Red = 90.0%, green = 99.0%, dark blue = 99.90% and light blue = 99.99%. 
Insets focus on energy ranges of interest. 
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Fig. 4. — Narrow Gaussian line (width < instrumental resolution) for a spectrum containing 
800 counts (left) and 1600 counts (right). Comparison of the detection limits, in equivalent 
width (keV), of the three methods over the energy band pass of Swift. The data are as 
follows; dotted blue - Bayes factor analysis, solid red - RS method, and dashed green - 
posterior predictive p-value analysis. 
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Fig. 5. — Broad Gaussian line (width = 0.2 keV) for a spectrum containing 800 counts (left) 
and 1600 counts (right). Comparison of the detection limits, in equivalent width (keV), of 
the three methods over the energy band pass of Swift. The data are as follows; dotted blue 
- Bayesian analysis, solid red - RS method, and dashed green - ppp . 
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Fig. 6. — Broad excess (width = 0.1 keV) for a spectrum containing 800 counts (left) and 
1600 counts (right). Comparison of the detection limits, in equivalent width (keV), of the 
three methods over the energy band pass of Swift. The data are as follows; dotted blue - 
Bayes factor analysis, solid red - RS method, and dashed green - ppp . Values below 0.7 keV 
have been excluded due to the width of the features being analyzed. 
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Table 1. Summary of spectral fits for all candidate spectra. 



Model" 



Photon 
index 



Line energy 
(keV) 



Line width 
(keV) 



Line norm. (xlO ^ 
photons cm"^ s~^) 



Equiv. width 
(eV) 



TVh 
(xlO^o 



cm 



GRB 050730 (T+692s to T+792s) 
1 
2 



2 03+° " 



1 98+°-^^ 

-'■•^°-0.13 



1.78 



+0.18 
-0.19 



r> 70+0.02 
'-'•''5-0.03 



1.14 



+0.48 
-0.44 



< Inst. res. 



0.34 



+0.08 
-0.16 



0.95 
22.5 



+0.58 
0.57 
+7.5 
-18.5 



50 
3400 



5.84^2-32 

800+6.57 
.»»_4 72 

+3.95 
4.03 



7.08 



64/55 
57/53 
47/52 



GRB 


060109 (T+109S 


to T+199s) 












1 


2.29l°-i? 










31, 


q+5.5 

•^-5.0 


2 


^•^■^-0.17 


74+0 03 

^■'^-0.03 


< Inst. res. 


1 1+3.08 
"^■-'^^-1.95 


94 


34, 


4+6.4 
■^-5.6 


3 


^•^'-'-0.18 


< 0.72 


'-'•^'5-0.06 


17.5t^3^8 


560 


39, 


•"-7.6 



48/42 

40/40 
35/39 



GRB 060111A (T+174S to T+234s) 
1 
2 

3 3.07t^;^^ 0.65 



o ncr+0.22 

O.UO_Q 20 

o.uy_Q 21 



64+° °^ 

U-"J*-0.03 



+0.09 
-0.06 



< Inst. res. 
< 0.13 



11.6 
14.0 



+15.2 

4.9 
+ 16.0 
-10.1 



73 
94 



29.91^5 
+5.7 
5.1 
+7.9 
-5.8 



32.6 
33.0 



52/50 
44/48 
44/47 



GRB 060111A (T+319S to T+339s) 
1 



1 07+014 

■■-•y ' -0.14 
1 04+014 



1.94 



+0.07 
-0.09 



0.79 
0.79 



0.02 
0.01 
+0.02 
0.01 



< Inst. res. 
< 0.15 



9.42 
9.42 



+4.86 

4.18 
+ 11.9 
-4.18 



80 
80 



18.811:^ 

19.5tt^ 



19.5 



+2.7 
-2.2 



69/61 
54/59 
54/58 



GRB 060115 (T+121S to T+253s) 
1 

2 1.85^ 



1.88l°i? 



1.82 



+0.12 
0.12 

+0.13 
0.11 



89+" "-'' 



0.81 



+0.07 
-0.07 



< Inst. 
0.10 



res. 

+0.06 
-0.05 



n fi7+0-44 



2.09 



+2.21 
1.22 



39 
100 



16.6t|? 

+3^1 
3.2 
+4.3 
-2.7 



16.3 
17.0 



93/80 

86/78 
82/77 



GRB 060124 (T-f 537s to T+542s) 

1 1.30+°;}! 

2 1.29tSlf 2.49^ 



1.13 



+0.19 

-0.24 



2.30 



+0.06 
0.01 
+0.21 
-0.23 



< Inst. res. 



0.48 



+0.17 
-0.11 



11.8 
57.7 



+6.3 
6.5 

+42.0 
-25.2 



800 

150 



29.5l|^ 



27 
18.3 



9.3 
7.8 
+9.7 
9.9 



72/47 
62/45 
51/44 



GRB 060202 (T-f 429s to T+529s) 

1 2.16i°;}J 

2 2.15tS|j 0.94^ 



2.12 



+0.10 
-0.12 



0.94 



1+0.03 
-0.02 
+0.05 
-0.08 



< Inst. res. 
< 0.34 



2.69 
4.94 



+ 1.51 
1.30 
+5.01 
-2.80 



54 

99 



47 

^' •■1-4.0 



48 



-1+4.8 
J^-4.2 



50.0 



+1.1 
-5.0 



109/103 
97/101 
96/100 
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Table 1 — Continued 



Model" Photon Line energy 


Line width 


Line norm. ( x 10 ^ 


Equiv. width 






index (keV) 


(keV) 


photons 


cm ^ s ^) 


(eV) 


(xlO^o cm-2) 




GRB 060210 (T+233S to T+353s) 

1 2 J2~^^'^^ 

2 2:71+°;! o.66tS:i 

3 2.68tVil 0.67tro! 


< Inst. res. 
06+° °^ 

U.UO_o.o3 




4.DO_2.4o 

7 nR+5-94 
7.O6I454 


63 
100 


20.6^1.9 
21.5+;^^ 
21.6tt° 


98/72 
85/70 
81/69 


GRB 060418 (T+119s to T+129s) 
1 1.82tr^l 

9 1 n9+0.13 9 /I9+0.02 
z i.oz_Q ^-^^-O.OS 
n -, 09+0.13 9 49+0.02 
1.0Z_Q ^•^^-0.04 


< Inst. res. 
< 0.14 




7 9-1 +3.00 
' •^-'^-2.97 
7 91+4.74 
' •^■'■-2.98 


190 
190 


24.3+_tl 
23.7tli 
23.7+t^ 


72/59 
56/57 
56/56 



GRB 060418 (T+169s to T+194s) 



1 


9 7n+0.22 
^- '^-0.19 










99 q+4.3 
^^■'-'-3.9 


62/52 


2 


9 (57+0.07 
-0.12 


^■"^^-0.02 


< Inst. res. 




58 


99 9+2.1 

z^.z_2.0 


52/50 


3 


-1 09+0.22 
i.oz_0 66 


< 0.75 


'-'•"3' -0.20 




2300 


ii.oiri 


43/49 



GRB 060428B (T+212s to T+252s) 
1 



O.UZ_Q -^6 

2 04+0- 18 

^•^^-0.16 



2.83 



+0.16 
-0.16 



0.77: 
0.76 



-0.03 
-0.02 
+0.05 
-0.06 



< Inst. res. 



0.09 



+0.05 
-0.03 



3.561J: 



8.38 



84 
87 
+4.92 
-3.69 



39 
100 



9.21 



+2.37 

-2.27 



78/63 
69/61 
63/60 



GRB 060428B (T+252s to T+418s) 
1 

2 2.58 



2 64+0- 14 

^•"^-0.14 



3+0.15 
-0.14 
9 00+0. 14 



6q+""2 

"•"^-0.03 

< 1.07 



< Inst. res. 

r, 00 + O.O6 



^•"^-0.37 
q 71+5.48 
■'■-8.19 



34 
2100 



2 08+1-63 
^•^^-1.50 
-1 70 + 1. 64 
'"^-1.51 



0.39 



+1.93 

-0.17 



58/64 
50/62 
48/61 



'^Models: [1] Absorbed power law, [2] absorbed power law plus a narrow Gaussian (width restricted to less 
than the instrumental resolution) and [3] absorbed power law plus a free- width Gaussian. Models containing 
blackbody components are not reported in this table as the fits were poorly constrained. All errors are 
quoted at 90.0% confidence. 
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Table 2. 


Summary of spectral fits 


for GRB 050822 (T+471 s 


to T+661 s). 






Model" 


Photon 


High energy 


Line energy 


Line width 


Line norm. 


Equiv. width 


TVh 






index 


cutoff 


(keV) 


(keV) 


(xlO^^ photons 


(eV) 


(xlO^" 








(keV) 














1 
2 
3 
4 


5.20t°il 

9 97+0.57 
^•^ ' -0.18 
9 9r+0.12 


L'-O-L-O.Ol 

n O9+0.03 
^•o^-o.oa 


< 0.47 






240 
300 


34.ot^:I 

< 12.2 

q K1+0.86 


90/40 
49/37 
70/39 
55/36 



'^Models: [1] Absorbed power law, [2] absorbed power law plus a free-width Gaussian, [3] absorbed cutoff 
power law, [4] absorbed cutoff power law plus a free-width Gaussian. 




2 4 6 8 10 



Energy (keV) 

Fig. 7. — GRB 050730 (T-f692s to T-f792s): RS method. Inset focuses on energy range of 
interest . 
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Confidence limits 
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Energy (l<eV) 

Fig. 8. — GRB 050730 (T+692s to T+792s): RS comparison between the absorbed power 
law models containing a single A^h component (wabs, dotted lines) and two components 
(wabs and zwabs, solid lines). Note that the feature becomes far less significant with the 
addition of the A^h column at the appropriate redshift. 
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Fig. 9.— GRB 060109: RS results. 
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Fig. 10.— GRB 060111A (T+164s to T+234s): RS results. 
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Fig. 11.— GRB 060111A (T+319s to T+339s): RS results. 
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Fig. 13.— GRB 060124 (T+537s to T+542s): RS results. 
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Fig. 14.— GRB 060202 (T+429s to T+529s): RS results. 




2 4 6 8 



Energy (keV) 

Fig. 15.— GRB 060210 (T+233s to T+353s): RS results. 
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Fig. 16.— GRB 060218 (T+159s to T+2770s): RS results. 
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Fig. 17.— GRB 060218 (T+2359s to T+2409s): RS results. 
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Fig. 18.— GRB 060418 (T+119s to T+129s): RS results. 
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Fig. 19.— GRB 060418 (T+169s to T+194s): RS results. 
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T ' 1 ' 1 ' r 
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Fig. 20.— GRB 060428B (T+212s to T+252s): RS results. 




Energy (keV) 

Fig. 21.— GRB 060428B (T+252s to T+418s): RS results. 
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Fig. 22. — Spectra modeled with an absorbed power law model: GRB 050730 (black), GRB 
060109 (red), GRB 060111A (T+174 s to T+234 s, green), 060111A (T+319 s to T+339 s, 
dark blue) and 060105 (light blue). Note the residuals around 0.7 keV. 
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Fig. 23. — Spectra modeled with an absorbed power law model: GRB 060202 (black), GRB 
060210 (red), GRB 060418 (T+169 s to T+194 s, green), 060428B (T+212 s to T+252 s, 
dark blue) and 060428B (T+252 s to T+418 s, light blue). Note the residuals around 0.7 
keV. 
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Fig. 24. — Spectra modeled with an absorbed power law model: GRB 060124 (black) and 
GRB 060210 (T+119 s to T+129 s, red). Note the residuals around 2.3 keV. 
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Facilities: Swift. 
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